#How Do Observers Assess Resolve?
#July 7, 2018
#BJPS Replication 3.R

#This file contains the replication code needed to replicate the elite benchmark analysis from the main text (located in the "Generalizing to Elites" section)

library(ggplot2)

### To generate the necessary dataframes, either run BJPS Replication 1.R, or run lines 10-15 below:
library(here) #Avoids needing to setwd()
dat3 <- get(load(file="Resolve Conjoint Cleaned No US 070718.RData"))

#Do parallel processing - speeds things up
library(snow)
cl <- makeCluster(8,"SOCK")

###

## Now, load elite benchmark data
dat.elite <- get(load(file="Elite Benchmarks.RData"))

#Regime type results from conjoint mass sample
conj2 <- subset(dat3, dat3$RegimeType!="Mixed")

set.seed(43215)
B <- 2500 #Number of bootstraps

### First, extract regime type results from each sample
#Regime type results from elite sample
knes.dat1 <- na.omit(as.data.frame(cbind(DV_A1=dat.elite$DV_A1, A=dat.elite$A))) #Drop missingness
boot.k1 <- rep(NA, B)
for (i in 1:B){
  j <- sample(1:nrow(knes.dat1), replace=TRUE)
  temp <- knes.dat1[j,]
  boot.k1[i] <- mean(temp$DV_A1[temp$A==1], na.rm=TRUE) - mean(temp$DV_A1[temp$A==0], na.rm=TRUE) 
}

clusterBootS2 <- function(data){
	i <- sample(unique(data$ID_numeric),length(unique(data$ID_numeric)),replace=TRUE)
	row.nums <- NULL
	for (j in 1:length(i)){
		row.nums <- c(row.nums, which(data$ID_numeric==i[j]))
	}
	return(data[row.nums,])
}

bootConjoint <- function(...){
	temp <- clusterBootS2(conj2)
	temp.res <- mean(temp$Chosen[which(temp$RegimeType=="Democracy")], na.rm=TRUE) - mean(temp$Chosen[which(temp$RegimeType=="Dictatorship")],na.rm=TRUE)
	return(temp.res)
}

clusterExport(cl, list("conj2", "bootConjoint", "clusterBootS2"))

system.time(conj.results <- parSapply(cl, rep(1,B), bootConjoint)) 
conj.results <- conj.results * 100

## Now, extract costly signalng results from each sample

set.seed(43215)

#Costly signal results from the elite sample
knes.dat2 <- na.omit(as.data.frame(cbind(DV_B1=dat.elite$outcmBW, Threat=dat.elite$Threat, Mobilize=dat.elite$Mobilize))) #Drop missingness
boot.k2 <- matrix(NA, B, ncol=2)
for (i in 1:B){
  j <- sample(1:nrow(knes.dat2), replace=TRUE)
  temp <- knes.dat2[j,]
  boot.k2[i,1] <- mean(temp$DV_B1[temp$Threat==1], na.rm=TRUE) 
  boot.k2[i,2] <- mean(temp$DV_B1[temp$Mobilize==1], na.rm=TRUE) 
}

#Costly signal results from conjoint mass sample

bootConjoint <- function(...){
	temp <- clusterBootS2(dat3) #Use dat3 instead of conj2 so we aren't dropping the mixed regime type observations
	temp.res <- c(mean(temp$Chosen[which(temp$curBehavior=="Public threat")], na.rm=TRUE) - mean(temp$Chosen[which(temp$curBehavior=="Nothing")],na.rm=TRUE), mean(temp$Chosen[which(temp$curBehavior=="Mobilized troops")], na.rm=TRUE) - mean(temp$Chosen[which(temp$curBehavior=="Nothing")],na.rm=TRUE))
	return(temp.res)
}

clusterExport(cl, list("dat3", "bootConjoint", "clusterBootS2"))

system.time(conj.results2 <- parSapply(cl, rep(1,B), bootConjoint))
conj.results2 <- conj.results2 * 100

## Now, combine dataframes
x_2 <- data.frame(y=c(boot.k1, boot.k2[,1], boot.k2[,2], conj.results, conj.results2[1,], conj.results2[2,]), z=rep(c("Leaders", "Public"), each=(B*3)), dv=rep(c(rep(c("Regime Type"),B), rep(c("Costly Signals"),2*B)),2), trt=rep(c(rep(c("Democracy"),B), rep(c("Public Threat"),B), rep(c("Mobilization"),B)),2))

#Collapse sample and treatment into a single variable
x_2$ztrt <- NA
x_2$ztrt[which(x_2$z=="Knesset" & x_2$trt=="Democracy")] <- "Knesset: Democracy"
x_2$ztrt[which(x_2$z=="US MTurk" & x_2$trt=="Democracy")] <- "US MTurk: Democracy"
x_2$ztrt[which(x_2$z=="Knesset" & x_2$trt=="Public Threat")] <- "Knesset: Public Threat"
x_2$ztrt[which(x_2$z=="US MTurk" & x_2$trt=="Public Threat")] <- "US MTurk: Public Threat"
x_2$ztrt[which(x_2$z=="Knesset" & x_2$trt=="Mobilization")] <- "Knesset: Mobilization"
x_2$ztrt[which(x_2$z=="US MTurk" & x_2$trt=="Mobilization")] <- "US MTurk: Mobilization"

##### Figure 5 #####
# Note: edit colors so it can print in greyscale
dev.new(height=3.5,width=10)
ggplot(x_2, aes(x=y)) + geom_density(aes(fill=z), alpha=0.75) + labs(x="Average Treatment Effect", y="Density") + scale_fill_brewer(guide=guide_legend(title="Sample"), palette="Set1") + facet_grid(~trt, scales="fixed") + theme_bw() + scale_y_continuous(breaks=NULL) + geom_vline(xintercept=0, linetype=3, alpha=0.5, show.legend=FALSE) 

#### Quantities of interest in text:
#Regime type:
round(mean(conj.results, na.rm=TRUE),digits=1) #Mass public: -4%
round(mean(boot.k1, na.rm=TRUE), digits=1) #Elites: -6%

#Costly signals:
round(apply(conj.results2, 1,mean, na.rm=TRUE),digits=1) #Mass public: 6.6% threats and 11.7% mobilization
round(apply(boot.k2,2,mean, na.rm=TRUE), digits=1) #Elites: 8.1% threats and 6.8% mobilization